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Abstract Weakly nonlinear analysis of a two dimensional sheared granular flow is 
carried out under the Lees-Edwards boundary condition. We derive the time dependent 
Ginzburg-Landau (TDGL) equation of a disturbance amplitude starting from a set of 
granular hydrodynamic equations and discuss the bifurcation of the steady amplitude 
in the hydrodynamic limit. 
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1 Introduction 

To control flows of granular particles is important in science and industry [T][2l|3]|3] . 
However, the properties of granular flow have not been well understood yet, because 
they behave as unusual fluids [5]. This unusual nature is mainly caused by inelastic 
collisions between granular particles. Indeed, there is no equilibrium state in granular 
materials because of inelastic collisions between grains, which suggests that granular 
materials are an appropriate target of nonequilibrium statistical mechanics [^. 

Although there are many studies of granular flows on inclined planes [71|^, the 
existence of gravity and the role of bottom boundary make the problem complicated. 
On the other hand, the granular flow under a plane shear is the simplest and an 
appropriate situation for theoretical analysis. Therefore, granular flows under a plane 
shear have been studied from many aspects such as the application of kinetic theory[21 
IIP) , shear band formation in moderate dense granular systems [lllll2j . long-time tail 
and long-range correlation function [T3l14 15 T6, 17, 18T9 , 20 , 2lll22). pattern formation 
of dense flow [23ll24|[251l26ll27ll28j , determination of constitutive equation for dense flow 
pgil3UU31| . as well as jamming transition [5^133 ll34ll551l5M?fj . 



Kuniyasu Saitoh 

Yukawa Institute for Theoretical Physics, Kyoto University, Sakyoku, Kyoto 606-8502, Japan 
Present address: Faculty of Engineering Technology, University of Twente, 7500 AE Enschede, 
Netherlands 

E-mail: k.saitoh@utwente.nl 



Hisao Hayakawa 

Yukawa Institute for Theoretical Physics, Kyoto University, Sakyoku, Kyoto 606-8502, Japan 
E-mail : hisao@y ukawa. kyot o- u . ac . j p 



2 



In this paper, we focus on the shear band formation in moderate dense granular 
gases observed in the discrete element method (DEM) simulations [11II12| . It is known 
that two shear bands are formed near the boundary and they collide to form one 
shear band in the center region under a physical boundary condition. A similar shear 
band formation is also observed under the Lees-Edwards boundary condition. Such 
a dynamic behavior of shear bands is reproduced by a simulation of granular hydro- 
dynamic equations [T5] derived from the kinetic theory for granular gases 38 S WiDl 
I41II42IIT3II44II45) . In addition, the linear stability analyses suggest that a homogeneous 
state of the sheared granular flow is almost always unstable [46II47II48II49II50I[5T1I52) . 

Amongst many papers, it is notable that Khain found the coexistence of a solid 
phase and a liquid phase of granular particles in his molecular dynamics simulation of 
a dense sheared granular flow [27II28| . He demonstrated the existence of a hysteresis 
loop of the difference of density between the boundary layer and the center region of 
the container by controlling the value of the restitution coefficient. It should be noted 
that the mechanism of an appearance of the subcritical bifurcation based on a set 
of hydrodynamic equations, differs from that observed in the jamming transition of 
frictional particles [53| . 

Recently, Shukla and Alam carried out a weakly nonlinear analysis of a plane 
sheared granular flow, where they derived the Stuart-Landau equation of a distur- 
bance amplitude under a physical boundary condition starting from a set of granular 
hydrodynamic equations |54II55|[5S] . They found the existence of subcritical bifurca- 
tions in both relatively dilute and dense systems, while the supercritical bifurcation 
appears in other parameter space. However, the Stuart-Landau equation cannot be 
used to explain the slow evolution of the spatial structure, because they adopted the 
method by Reynolds and Potter [57] which does not include any spatial degrees of 
freedom. We also indicate that their perturbation is based on the analysis for a finite 
size system, in which the relation between the perturbation parameter and shear rate 
becomes unclear, because the shear rate is fixed to unity in their paper. 

In this paper, we derive the time dependent Ginzburg-Landau (TDGL) equation 
under the Lees-Edwards boundary condition ^58, as a spatially dependent amplitude 
equation of the disturbance fields starting from a set of granular hydrodynamic equa- 
tions |59ll60ll6l1l62ll63ll64| . To reduce the number of control parameters, we only focus 
on the behavior in the hydrodynamic limit. We discuss the bifurcation in the hydrody- 
namic limit from the results of the coefficients of the TDGL equation. The organization 
of this paper is as follows. In the next section, we explain our setup and basic equations 
of a two dimensional sheared granular flow. In SecO we summarize the results of the 
linear stability analysis. Section [4] is the main part of this paper, in which we derive 
the TDGL equation with the aid of the weakly nonlinear analysis. Finally, we discuss 
our analysis and describe our conclusion in Sec[5] 

2 Setup and basic equations 

Let us introduce our setup and basic equations. To avoid difficulties caused by physical 
boundary conditions, we adopt the Lees-Edwards boundary condition, in which the 
upper and lower image cells move to the opposite direction with the speed 17/2 [58) . The 
geometry of our setup is illustrated in Fig[l]with the Cartesian coordinate x = {x,y). 
Because we adopt the diameter of a granular disk d and U/2 for the unit of length and 
speed, respectively, the shear rate U/L is reduced to e = 2d/L in this dimensionless 
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Fig. 1 Geometrical setup of a two dimensional sheared granular flow under the Lees- Edwards 
boundary condition. The upper and lower image cells move to the opposite direction with the 
dimensionless speed 1. The dimensionless width and height of each cell are W/d and L/d, 
respectively. 



unit. In the following, we also use the mass of a granular disk m and 2d/U as the unit 
of mass and time, respectively. 

We employ a set of hydrodynamic equations derived from the kinetic theory of 
granular gases [44]. Although the angular momentum and the spin temperature are 
included in the hydrodynamic equations, we ignore such rotational degrees of freedom 
to simplify our analysis. If the friction constant is small, this simplification can be 
justified, because the effect of the rotation of granular particles during the collision can 
be absorbed in the normal restitution coefficient [651166] . 

We present the derivation of the following set of dimensionless hydrodynamic equa- 
tions in [Appendix A| 

(9t + V • V) = ~uV ■ V (1) 

1/ (a* + V • V) V = -V • P (2) 

(!//2)(at-f V- V)6» = -P : Vv- V-q-x , (3) 

where v,w = (it, w), 6, t and V — {d/dx, d/dy) are the area fraction, the dimensionless 
velocity fields, the dimensionless granular temperature, the dimensionless time and the 
dimensionless gradient, respectively. The pressure tensor P, the heat fiux q and the 
energy dissipation rate x a^re given by 

p = [p*(u)e - C{y)0^'^ (V ■ v)] 5,, - i^*(v)e^'^D[^ , (4) 

q = -K* {u)e^^^Vd - X* (y)&^l'^^u , (5) 
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1 - 2 / Nai/2 
4v27r 



46»-3^/|e^/^ (V- 



(6) 



respectively. Here, D[j {i,j = x,y) is the deviatoric part of the strain rate 

D',j = i {\/jv, + \/,Vj - S,j\/ ■ v) , (7) 

and p*{u)e, i*{v)e'^/'^, ri*{u)6^^^, k*{u)6^^'^ and X*(ty)e^^'^ are the static pressure, the 
bulk viscosity, the shear viscosity, the heat conductivity and the coefficient associated 
with the gradient of density, respectively. The explicit forms of them are listed in Table 
[TJ where we adopt 

{l-v) 
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Table 1 The functions in Eqs.l(4j-(|6ll, where e is the restitution coefficient. 





lu[l + {l + e)ug{u)] 




CM = 






ri*{u) = 


r^\a(.i^)'^ 1 (l + e)(3e + l),, 1 / (l + e)(3e-l) , 1 A 
V 2 [ 7-3e ' 4(7-3e) ' V 8(7~3e) tv J '^-^ 


f e)!/23(,,)j 


K*{u) = 


./T^-r 9(^)-l 1 3(2e^+e + l) / 9(l + e) (2e~ 1) 
^'^^ [(l + e){19-15e) ' 8(19-15e) ' \ 32(19-15e) ~^ 




A*H = 







for the radial distribution function at contact which is only valid for u < 0.7 [6711681 
IW, 70 . It should be noted that the expression of niu) in Ref. |12) contains an error (see 
[Appendix A[ ). 



3 Linear stability analysis 

In this section, we present the linear stability analysis of a sheared granular flow under 
the Lees-Edwards boundary condition. Although the analysis is essentially same as 
those in the previous studies 46.47.48.49.,50.,51.i5^ . it is necessary as the basis of the 
weakly nonlinear analysis. 



3.1 Linearized equation 

We introduce the hydrodynamic field and the homogeneous solution of Eqs.([l])-([3l) 
as (j> = {i',u,w,9)'^ and <^o = {i^Ot ^V^O, 9q)'^ , respectively, where the upperscript T 
represents the transposition, uq is the mean area fraction and 



2 (1 - e2) u^g{uo) 



(9) 



is the mean granular temperature. Thus, in the hydrodynamic limit e ^ 1, 1 — is 
scaled as 1 — = with the fixed ^o- The disturbance field is defined as 4>{x,y,t) = 
(p — (po which is transformed into the Fourier series 

kyQ 7^0 kyQ 

where the upperscripts L and NL respectively represent the layering mode {kx ~ 0) 
and non-layering mode (kx 7^ 0), and with I = L or NL is the amplitude. The 
so-called Kelvin mode is defined as 

k(f) = {kx,ky(t)) = {kx, kyo - etkx) , (11) 

where kyQ = ky{0) and the coefficient 4'\^(^t) defined with the imaginary unit i as 

't'i.it) = ('^k(t).«^*k(i).*™k(i).^k(t))''' ■ (12) 
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We also introduce <^k(t) = (^k(t)i '^k(t)i ''^k{t)j ^k(t))'^ f"'" the convenience of the analy- 
sis. If we linearize Eqs.lT}-©, 'fi]s.{t) satisfies 



dt 



(13) 



where the convective term is canceled because of the Kelvin mode Eg. pip . The time 
dependent matrix C{t) is decomposed as 



= jCo{kx, kyo) + tjCi{kx, kyo) + t C2{kx, kyo) ■ 
The matrices Co{kx, kyo), Ci{kx, kyo), C2{kx, kyo) are respectively given by 



(14) 



i^oikx, kyo) = 

/ vokx vakyo \ 

^"^^yO ~ Pq^x ~^okx ~ ^k —^okxkyQ — e ^"^^kyQ — ^kx 
efkx - p'okyo -^okxkyo -Cofe'o - f fc' ^^kx - ^kyO 

V e^ci - 2Aofc^ C2kx - 2erjQkyo C2kyo - 2erjokx e^c^ - 2^0^^ / 

/ -euokx \ 



C-l{kx, kyo) = 



£-2{kx, kyo) 



e^^kx eriokxkyo eCo^x 



'^kx 



ep'okx e^okl e{2^o + rio)kxkyo ^^kx 
\4:eXQkxkyo 2e^riokx —eC2kx Aenokxkyo/ 

/ ° \ 

-e^^fc^ 

-e2(Co + f)fc' 

) -2e^KoklJ 



\-2e^\ok^ 



(15) 



(16) 



(17) 



where k = kx + fc^Q and ci , C2 and C3 are respectively given by 



/2/ , N/)3/2 
Cl = ??l - \/ -(30 + !^05l)Po ' 

V TT 

C2 = 2po - -e I'oSO^O , 



C3 



??0 
200 



1/2 



(18) 
(19) 
(20) 



The explicit forms of the coefficients of the Taylor expansion, i.e. gQ,PQ, (,q, r]Q, kq, Xo,gi 
and Po are respectively given by Eqs.([Sll)-((Hni), ((M]) and (|100p in [Appendix B| 



3.2 Non-layering mode 



The solution of Ea. p3p is obtained by the parallel procedure in Refs. [19ll2UII21) for 

the case of the non-layering mode (kx 7^ 0). In [Appendix C| we perturbatively solve 
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Ea. (|13|l by scaling the wave number as k(t) = eq(t) and find the components of '/'q(t) 



"q(t) 



w, 



q(t) 



where we defined 



et 



1 



-j=l==E^'\t)sinu;{t) 
rf „(3), 



9^^t) = ^jE^'Ht)+^-fE^'\t) cos.it) 



E^'^'it) smuj{t) 



(21) 
(22) 

(23) 
(24) 



and the frequency 



E^'^\t) = exp 



£^^^(t) = exp 



E'^^\t) = exp 



4 

2 / 2, e ,3 

-'7j;'-6 ( e * + 

4 

2 / 2, e ,3 



etdl + (et)2 + In <^ + i/l + (et) 



(25) 
(26) 
(27) 

(28) 



The positive constants J, ra, and rc are respectively given by Eqs.([T22]) and (|123|) 
in [Appendix C| 

From Eqs.(l2lJ-((24ll, tp^^-^ decays to zero in the long time limit as indicated in the 
previous works 48,50 . Therefore, the nonlayering mode is linearly stable. It should be, 
however, noted that H^^t) involving the convective effect is only necessary for 7^ 
[19II20II21) . Thus, we can solve Ea. (|13|l for g^ = separately in the next subsection. 



3.3 Layering mode 

In the case of the layering mode (kx = 0), Ea. (|13|) is reduced to the eigenvalue problem 

£0(0, fcj^o)'/'fc^o ='^{kyQ)ifi\^ , (29) 

where a{kyQ) and '■p\^^ are the eigenvalue and eigenvector of £(0, fcyo), respectively. 
We also define the left eigenvector (p\ ^ as 



'^fc„o'^o(0, feyo) = cr(fcyo) 



(30) 



In [Appendix D| we perturbatively solve Eqs. (|29|) and (|30|) with the scaling kyo — eq 
and find the dispersion relation 



'^(9) = £^(''29^ + r4,q^) 



(31) 
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q 

Fig. 2 The dispersion relation cr{q)/e'^ , where the open circles and the solid line represent the 
numerical results and Eq. |I31II . respectively. The maximum value is given by the scaled wave 
number qc = 0.057. Here, we used e = 0.01, e = 0.99 and vq = 0.4. 



which is maximum at qc = \/— r2/2r4, where r2 and r4 are given by Eas. ()159|) and 
(|161|l in [Appendix D| The right and left eigenvectors are respectively given by 

/ (1) ' 

ipq = [iyq,Uq,Wq,eq) =1 ^-J- , a\ \ eC^^ , —j , (32) 



, ,-(1) „ -(1) -(1) ,(1) 

(fig = {Uq,Uq,Wq,eq) = | —jj^ j , ',eC^^, ^^^^^ H — ] , (33) 

where a^^^ = (aJ^\a^^V. a^^^ = (ai^^4^^)'^, C'-Pi and are given by Eas. (fT57)) . 
(fT68ll . (IT63ll and (flTOll in [Appendix D[ respectively. 

Figure [5] shows the dispersion relation a(q)/e^ , where the open circles and the solid 
line represent the numerical results and Eq.(|3T}, respectively. In numerical calculation, 
we solved the eigenvalue problem Eq.^ by LAPACK [TT] with e = 0.01, e = 0.99 
and z/Q = 0.4. In this figure, the maximum value 

Oc = = r2qc + r^qc 34) 

is given by qc — 0.057. It should be noted that the imaginary part of cr(g) is always 
zero. 



4 Weakly nonlinear analysis 

The linear stability analysis is only useful to know whether the considered base state 
is stable. If we are interested in the structure formation after the base state becomes 
unstable, we need, at least, a weakly nonlinear analysis. Let us introduce a long time 
scale and long length scales as r = f^t and z = (^, C) = ^{x,y), respectively, to 
characterize the slow and large scale evolutions of structure. Thus, the derivatives are 
replaced by 

dt = e^dr, V = e(a^,9^). (35) 
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The slow evolution of hydrodynamics variables are obtained from the evolution 
of the neutral solution of the linearized equation. The neutral solution at the most 
unstable mode qc = (0, qc) is given by 

0„ = ^L(C,r)0^^e''«=';+c.c. , (36) 

where each component of is the corresponding one in Ea. p2|l aX q — and 
c.c. represents the complex conjugate. It is notable that the amplitude A^{C,,r) is 
independent of ^, because the non-layering mode qx ^ Q are linearly stable. Thus, if 
we adopt the conventional approach in which the amplitude equation is obtained from 
the expansion around the neutral solution, we cannot discuss the structure evolution 
in ^ direction. 

If we carry out the weakly nonlinear analysis using (^n, the amplitude equation 
for A^{C,,t) only depends on but the disturbance in the ^-direction also exists in 
the two-dimensional granular shear flow. Let us try to introduce a hybrid approach to 
involve ^ dependence in shear flow. For this purpose, we may rewrite (^(a:, y, t) in Eq. 
P0[) in the vicinity of q = qc as 

<^„~a(e,C,r),^^,e''l(")-" + c.c. , (37) 

where the wave number q(r) involve the contribution of the deviation qc, i.e. q(r) = 
qc + 5(\{t). In addition, we need to include the contribution of the non-layering mode 
'^qM ~ ('^qCr)! *'^q(T)) *™q(T)i ^q(T))'^ whcu wc are interested in the case of qx / 0. 
Thus, Ea. (|37l) may be replaced by the hybrid solution 

k = C, r)c^\. + ^""^C, C, T^ll^^y^^^-^ + c.c. 

- ^(?, C, r)[<^^, + <^^(';)]e^'i(")-^ + c.c. , (38) 

where we have used a strong assumption that the amplitudes a(^, r) and yl^'^(^, r) 
are scaled by the common amplitude A(^,^,r). If we carry out the weakly nonlinear 
analysis using (j>\^ instead of 4>n, the TDGL equation might depend on ^. Strictly speak- 
ing, we cannot justify the above hybrid approach between two different modes, i.e., the 
layering mode and the non-layering mode. Nevertheless, we will take into account ^ 
dependence in the TDGL equation phenomenologically. 

Now, let us proceed the explicit calculation of weakly nonlinear analysis. To avoid 
the confusion from the uncertain part in the hybrid approach, we first derive the one- 
dimensional TDGL equation in Sec l4.1l for only the layering mode, and give the hybrid 
TDGL equation in Sec l4.2l including the contribution from the nonlayering mode. 



4.1 Weakly nonlinear analysis of the layering mode 

In this subsection, we derive TDGL equation as the amplitude equation for the layering 
mode at the most unstable wave number. This subsection consists of three parts. In the 
first part, we expand the amplitude A}^{C,,t) and the matrix £o(0,egc) introduced in 
Ea. H15|) . In the second part, we will derive TDGL equation at 0{e^) which is sufficient if 
the bifurcation is supercritical. In the last part, we will present higher order calculation 
which is necessary if the bifurcation is subcritical. 
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4.1.1 Expansions of amplitude and matrix 



In this part, we prepare the expansions of the amplitude and the matrix in terms 
of e, which is necessary for the weakly nonlinear analysis. From the straightforward 
calculation, A^{(^,t) and Co{0,f.qc) can be expanded as 



£0(0, egc) = eMi + e^M2 + ... , 
where the matrix jCo is introduced in Ea. (|15|l and 



Ml 



( 



/ v^qc \ 
0-1 

-P'oqc -^qc 

\ 2pogc / 







M2 = 



2 1c 2 1c 



-(Co + f )9c' 









C3 - 2nQql) 



(39) 
(40) 



(41) 



(42) 





\ ci -2r;ogc 

Substituting Eqs. (|36P and H39p into Eqs.((T])-((31) and collecting the order of e, we can 
obtain a series of terms of equations. 

Multiplying the left zero-eigenvector of £o(0, e^c), we will obtain the amplitude 
equation. It is not easy to obtain the left zero-eigenvector in general, but fortunately 
<^g^ introduced in Ea. H33|) plays a role of the zero-eigenvector in the limit e — >■ thanks 
to Eqs.([29]) and 



4.1.2 The TDGL equation at O(e^) 

The first nonzero terms appear at O(e^), where the coefficient of e^'^"'' satisfies 

where </3g_, is introduced in Eq.((32|). At 0{e^), the coefficient of e*'^'' satisfies 

^^JrA\ = M2^^Ai+'^d^c^\ +^3A\\A\f , 
where D and As are given by 



V = 





(Co + %/2)Wg, 



\ 



( 



-P2^t - (Pl +P2)vljqJeQ 
\ 2Vq^Wq^{pieqJeQ+p2Vq^) ) 



(43) 



(44) 



(45) 



respectively. 

If we multiply the left zero-eigenvector to Eq. (|44p introduced in Ea. H33p . we 
obtain the TDGL equation: 



(46) 
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where we have used the normahzed condition (p]^^Lf^^ — 1, and d and /3 are given by 

Tig ~ 

d = —{uq^Uq^ + Wq^Wq^) + 5o™«c™«c + "^l^odqaSqa > (47) 
P = 2eq^Uq^Wq^ [p^^q^ ^ ^'^"^ ^ {^'^"^^ ^ '"'^ ' ^^^^ 

respectively. 

Substituting Eqs.dSl]) and (gHJ to Eqs.lgTI and (gSJ, the leading terms of e give 

d = J , /3 = e/3 , (49) 

where d and /? are listed in Tabled It is notable that the coefficient /3 becomes higher 
order of e. Therefore, we need to rescale the amplitude as 

A^,{(;,r)^e'/^A\iC,r) (50) 

and the TDGL equation for A^i((,t) is reduced to 

drA^i = acA^i + dd'^A^i + pA^i \A^i f . (51) 

The scaling relation Ea. (|50|l indicates that the amplitude of is extended as 

e'/^A^l (C, r) + e^/^A^2 (C, r) + e^^^A^3 {C,r) + ... , (52) 

where A^j = e^^^A^ (j = 2, 3, . . . )• Thus, converges to zero in the hmit e -5> 0. 

Let us compare Eq. (|49p with the numerical result, where we solve Ea. (|29p by LA- 
PACK and calculate /3 from Eq. (|48p . We find ctc and d are always positive and Ea.(|49[) 
perfectly agrees with the numerical results (FigO. We find /? < in < < 0.245, 
thus, a supercritical bifurcation can be observed in the dilute regime. On the other 
hand, and a subcritical bifurcation , i.e. > appears in > 0.245. 

It should be noted that there is no hysteresis behavior even for the subcritical 
bifurcation. Figure |4] is a schematic image of the subcritical bifurcation of l^'^il, where 
a hysteresis loop is realized by the paths (i), (ii) and (iii). Because we restrict our 
interest to the case of e > from the definition, the paths (i) and (iii) cannot exist. 
Therefore, such a hysteresis behavior cannot be observed in the hydrodynamic limit. 

4.1.3 Higher order expansions 

Because of ^ > in i^o > 0.245, we need to proceed our calculation to the higher order 
expansions. At 0{e'^) and O(e^), the coefficients of e^'^"'' satisfy 



^q.drA^ ^M2^iA!^ + m^^A!^+J^3{Ai A!^ +2A!^\A'^f), (53) 
V^.drA^ = M2'fi^^Al + VdlA\ + mAYAf + 2A\\a\\^ + aY" aY + 2\a\\^a\) 



+Af5A'^\A\:\'' + B{A\: d'^A\: +2\A\:\'d'^A'^)+c{A\: {d^A\:y + 2A\:\d(;A'^n, (54) 

respectively, where A^ {j = 1, 2, 3) represents the complex conjugate of A^ and 



~2p'^uleqjeo 

Wal^qJqaWqa/Oo/ 



(55) 
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Although the vectors B and C can be written exphcitly, we do not need these analytic 
forms in later discussion. 

Let us introduce the envelope function 

A^{C,T) = A\{C,T) + eA^{(;,r)+e^Al{(:,T) , (56) 

which is used by many authors to derive higher order amplitude equations [7211731171] . 
Summing up Eas. (|44)) . ((53)) and ((54)) . we obtain 

\B{A^'^dlA^* + 2\A^fdlA^) + C{A^\dc^A^f + 2A^\dc^A^\^}^ . (57) 
Then, multiplying (p]^^ to Eq.((57l) we find 

drA^ = (TcA^ + dd^A^ + l3A^\A^f + eSAL|i'L|'* 

[fe(i'L^92^~L* _^ 2\A^fd'^A^) + c{i'L*(a^i'L)2 ^ 2i'L|9^i'L|2}j , (58) 

where b = (p\^B, c = <^g^C and 

^ = ■^'^qJlci'^PsdqcWq, - PZ^q.Vq,) ■ (59) 

Substituting Eqs.(I32I and P5)l to Eg. (15^ . the leading terms of e give 

7 = £7 , 6 = , c — e'c , (60) 

where 7 is given in Table [J] Although 6 and c can be written explicitly, we do not 
need such analytic forms in later discussion. It is notable that the coefficient 7 becomes 
higher order of e by substituting Eas. H32p and (|33p . Thus, the rescaled envelope function 
^^{C,t) = e^/^A^{C,T) satisfies 

drA^ = acA^ + ddlA^ + ^iL|^L|2 e^Ah\A^f 



where the TDGL equation including the term of A'^IA^I'* is given in the first line. 

Let us compare Ea. (|6Up with the numerical result, where we solve Eq. (|29[) by LA- 
PACK and calculate 7 from Ea. (|59p . Figure [3] exhibits a complete agreement between 
Eas. ((59)l and ((60)) . From this result, for 0.245 < vq < 0.275, the growth of disturbance 
is inhibited by the nonlinear term e7^'^l^'^l^ finite steady amplitude can be ob- 
served. For > 0.275, we need to calculate higher order expansions, however, it is too 
complicated to perform in this paper. 
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0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.1 0.2 0.3 



t'O T^O 

Fig. 3 The TDGL coefficients /3 and 7 for (a) < < 0.7 and (b) the dilute regime 
< < 0.3, where the open circles and the open squares represent the numerical results of 
/3 and 7, respectively. The solid and the broken lines represent the analytic results of jB and 7, 
respectively. Here, we used e = 0.01. 




e 



Fig. 4 A schematic image of the subcritical bifurcation of where a hysteresis loop is 

realized by the paths (i), (ii) and (iii). 



4.2 Hybrid approach of weakly nonlinear analysis 

In the previous subsection, we have obtained the amplitude equation for the layering 
mode. The derivation is straightforward and the obtained amplitude equation has a 
reasonable form. The equation, however, only depends on ^, and thus, we cannot dis- 
cuss the spatial structure along the mean flow direction ^. To improve this unsatisfied 
situation, we adopt the hybrid approach as mentioned, though it is hard to justify this 
approach. Fortunately the contribution of f^jl,-^ except for the diffusive mode becomes 
irrelevant as time goes on. Thus, (^^^ still can play a role of the left zero-eigenvector in 
our calculation. 

Let us expand the amplitude of into the series of e as 

A{^,C,r) = eAi{^,C,r) + e^A2{^,C,r) + e''A3{^,C,r) + ... . (62) 

If we use instead of tpn and multiplying the approximate left zero-eigenvector i^^^ , 
we obtain the TDGL equation of j4i(^,(",t) as 



drAi^acAi+di{T)dlAi+d2{T)d^d^Ai + ddlAi+l3Ai\Ai\^ , (63) 
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Table 2 The cxpHcit forms of d, /3, 7, di{T) and 1^2 (t). 
d = ^aWaW + §^ (pod['^ + uoeoJg.-4'^) 

^ = -PoPi) (pofli^' + i^o^o Jgc4^^) + ■^'3cC'^iPo{poP2 - (pi +P2)Po}4^']4^' 

Mr) = {(f +«o)4''«q(r) + e;^(po4''+'^o9oJgc4'')eq(,)} 

/{l + j^(Po4^'' - 2poJgc4^^)'^q(T) +4^'"q{T) + BoJ'^qa + ^'ofo Jgc4^' )f q(r) } 

'?2(-r) = ?o4^'«'q(T)/ {1 + j5^(Po4^' - 2poJlJc4^')'^q(T) +a^^'wq(r) + 9j^7I^(Po4^' + l^O^O ^qc4^' )^'q(r) } 



where we introduced the time dependency in the diffusion constants 

d-lir) = {Y(i'9c^tq(r) +*gc«^q(T)) +Coitgc^q(T) +2«:o%^'q(T)} 

/{l + ^l^^l^)), (64) 

rf2(T) = Co(w<;o"'q(r) +*<?c'"q(r))/(l +'^9o'^q(r))> (65) 

which decay to zero because of Eqs. (|21|) - H24|) . To obtain H63p we ignore contributions 
from Vq(ii-) except for the diffusion coefficient, because they exponentially decay to 
zero. 

Substituting Eqs. (|21I) - (I24I) . H32p and H33p . and the introduction of the scaled am- 
plitude 

Ai(C,C,r)EEel/2^i(C,C,r) , (66) 
we find the TDGL equation of Ai(^, r): 

drAi = acAi + di{T)dlAi + d2iT)d^di;Ai + dd^Ai + pAilAif , (67) 

where the explicit forms of di (r) and d2 (r) are listed in Table [J] 

From the parallel argument to obtaining Ec ljOf |) . the scaled envelope function 
introduced 

A{^, C, r) = 6l/2{Ai(C, C, r) + eA2i^, C, r) + e^Asi^, C, r)} , (68) 

satisfies the TDGL equation of y4(^,(^,r): 

drA^acA + di{T)dlA + d2{r)d^d(;A + ddlA + l3A\Af + ejA\Af , (69) 

where we truncated the higher order terms of 0(e'^) in Ea. H69p . 

Figure [5] shows the time evolution of di(T) and d2{r), where the analytic results 
perfectly agree with the numerical calculation of Eqs. (|64I) and (|65p based on LAPACK. 
Because diij) and d2(r) decay to zero, Eqs. (|67P and (I69p respectively reduce to Eq. (|51P 
and the first line of Eq.((6T]) in the long time limit. This result is consistent with the 
observation in the simulation 12 in which the shear band finally becomes parallel to 
mean-flow direction, though the mathematical justification of our hybrid approach is 
difficult. 
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Fig. 5 The time development of di{T) and <i2(''-). The open circles and open squares represent 
the numerical results of diir) and d2{T), respectively. The solid and broken lines respectively 
represent the analytic results of di(T) and d2{T). Here, we used e = 0.01, e = 0.99 and uq = 0.4. 



5 Discussion and Conclusion 

Let us compare our results with the previous studies |54II55I[56] . The previous studies 
only derived Stuart-Landau equation which is independent of the position, while we 
obtain TDGL equation which can discuss the slow evolution of long-wave spatial struc- 
ture. We have demonstrated that the coefficient of j4|A|^ can be calculated explicitly 
based on a systematic perturbation method in terms of small e, which has not been 
achieved by previous studies. The appearance condition of the subcritical bifurcation 
is slightly different from that of the previous studies. We believe, however, that the 
result becomes similar to that of the previous studies, if we analyze a finite size sys- 
tem around most unstable mode. On the other hand, it is hard to justify our hybrid 
approach to introduce the time dependent diffusion coefficients di{r) and d2(T) in the 
TDGL equations Eas. (|67p and (|69p . though the result seems to be reasonable. The 
mathematical justification of the hybrid approach will be our future work. 

In conclusion, we have derived the TDGL equation starting from a set of granular 
hydrodynamic equations. From our results, we find the homogeneous state is always 
unstable and a supercritical bifurcation can be observed in the dilute regime < i^o < 
0.245. On the other hand, a subcritical bifurcation is predicted in vq > 0.245 and we 
find the amplitude of disturbance can be converged by the nonlinear term e7A|A|^ in 
the range 0.245 < < 0.275. In the case of i^q > 0.275, higher order expansions are 
necessary, however, such calculations should be performed in a future work. 
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Appendix A Derivation of the coefflcients in Table 111 

In this Appendix, we derive the coefficients in Tabled by using the dimensionless quantities 
based on the kinetic theory 1441 . At first, the energy sources Xaa and X3 for smooth disks 
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-1 in Ref.gi]) are 



2dP 



\%T -Z-K^I'^dT^I'^{V -u) 



(70) 



and X33 = Oi where d, T = m((c — u)'^)/2 and u = Uv/2 are the diameter of a disk, the 
granular temperature and the velocity field, respectively. It should be noted that we adopt the 
different definition for the granular temperature T from Ref. |44l to keep the dimension of the 
energy. In Eg. HTOI I. the bulk viscosity 5 is given by 



4m 



(l + e)i.2gHTi/2 



(71) 



where u and g{u) is the area fraction and the radial distribution function at contact, respec- 
tively. Thus, the factor of Eq. jTOj is given by 



9MT 



1/2 



(72) 



If we introduce the mass density of the system p = Amu/ {nd'^) and the mass density of a disk 
pp = Am/ (ird^), Eg.l lTOII is reduced to 



2pp7rl/2d' 
and the energy loss rate is given by 



29(i.)Tl/2 [sT - 37rl/2dTl/2(V ■ u)] , 



X ■ 



Xaa _ 1 - 



2g(i.)Tl/2 8T - 37ri/2dTi/2(V ■ u) 



The pressure tensor is given by Pij = pTSij + paij + Oij , where 

2mTi/2 

"''^^ = - .V2,,(.)(5-3.) + - ' 

Oij = (2pTug{u)r - • \i)5ij - ^D-^ + i/g(u)rpaij 

with r = (1 + e)/2. Then, wc find 

Pij = [P - CV • u] <5ij - + [1 + r!^g(i/)] pa,;j , 



(73) 
(74) 

(75) 
(76) 

(77) 



where the static pressure is given by p = pT [1 + (1 + e)i/g{i')]. The second and third terms on 
the right-hand-side of Ea.l l77ll can be rewritten as ^D'^^ — [1 -I- ri'g{i')] paij = V^'ip where the 
shear viscosity r) is given by 



4mTl/2 
7rl/2d 



g(v)-i , (l + e)(3e+l) 



7- 3e 



4(7 - 3e) 



(l + e)(3e- 1) 1 



8(7 - 3e) 



(l + e)v^g(v) 



(78) 



Therefore, we find Pij = [p - ■ u] &ij - r]Dlj. It should be noted that Eq.(70) in Ref. [ill 
should be multiplied by rm. The translational energy flux is given by qa = padjjj /2 + 0app /2, 
where paapp/2 and Octpp/2 are given by Eq.(89) and (100) in Ref. |44| . respectively. From 
Eq.(lOO) in Ref. [44], we rewrite qa as 



1 3 1 

qa = -paapp - ?VT + -rug{u) ■ -paapp 



1 + -rug{u) -paapn - ^VT , 



We introduce Kp and Ap as ^pa^pfj = —KpWT — ApVp, where 



4mTi/2 



ag{v)r{n - 15r)7ri/2 



1 + -ug{u)r'^{Ar-S) 



^a^^/^{2r-l){l-r) ^.,/^ d(u^g{u)) 



2vg{v){n - 15r) 



dv 



(79) 

(80) 
(81) 
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If we write the energy flux Qc as Qa = — kVT — AVp, we obtain the heat conductivity k 

/ 3 , A , 16mTl/2 r gM'^ 3(26^ + e + 1) 

K = K„ 1 + -rvg(v) ] + i= h -^^ -u 

'^V 2 ^^^7^^ o-7rl/2 L(l + e)(19- 15e) 8(19 - 15e) 

9(l + e)2(2e- 1) 1 1 , ^ 9 , si 



32(19 - 15e) 47r 
and the coefficient associated with the gradient of density A 

f 3 , A 3o-7ri/2e(l - e) r . ^ 1 , si ^d{v'^g(v)) , , 

A = Ap 1 + -rug{u) = \ ' [Ag(,u)-^ + 3(1 + e)u\ - ^ '> T^'^ . (83) 

V 2 / 8(19 — 15e) u du 

We should note that the third term on the right hand side of Eq. |I82II differs from our paper 
1121 . Indeed, the coefficient l/47r in the last term on the right hand side of Eq. II82I I is different 
from l/27r. 

Now, we non-dimensionalize the static pressure, transport coefficients and the coefficient 
associated with the gradient of density with the aid of m, d and [7/2 as 

K = ppdUK*(u)e^^'^, A = d^^j x*{u)e^^^, 

where p*{u), ti*(u), k,*(u) and X*{u) are dimensionless quantities listed in Tabled 



Appendix B The Taylor expansion of the functions in Table [T] 

The functions in Table [T] are expanded into the Taylor series as 

gM = go + gii^ + 921^^ + ■ ■ ■ , (84) 

eoi>~^p*{i>) = PQ+piu + P2u'^ + ... , (85) 

ey^'H'iy) = «o + iiv + ii2y^ + . . . , (86) 

ey^v'^m^u) = r]o+ riiiy + r)2U^ + ... , (87) 

dQ^^I/~^K''{h>) = Kq + KlU + K2l^^ + ■ ■ ■ , (88) 

eQ^^u-^\*{u) = Ao + Xiu + X2U^ + ... . (89) 
Similarly, the derivatives are also expanded into the Taylor series as 

dv 

9^--^^^ = + + C^.^ + ... , (91) 

el'\-^'^=v', + v{u + ^'2U^ + ... , (92) 
dp 

el^^u-^ '^'^J''^^ = k[, + k[u + k'2u'^ + . . . . (93) 

In the following, we show the explicit expressions of the coefficients which are used in the text. 
The coefficients associated with the radial distribution function are given by 

_ 25 - Tup _ 34 - 7uo _ 43 - lup _ 52 - Yup 

16(1 - ^p)3 ' 16(1 - up)^ ' 16(1 - 1/0)5 ' 16(1 - vo)6 • ^ ' 

The coefficients associated with the static pressure are given by 
Pl = ^(1 + e.)igo + fogi)Sp, P2 = ^(1 + e)(3i + upg2)dp, ps = i(l + e){g2 + ypg3)dp . (95) 
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The coefficients associated with viscosity are given by 



where 

2V 8(7 -3e) ' TT J ' ' V2 7-3e 

The coefficients associated with the heat conductivity are given by 



where we have introduced 



1 



(96) 



.^(l+e)(3e-l)^lX ^^^^^ /| 1 _ ^^^^ 



V 32(19 -15e) 47ry ^ ^ ' (1 + e)(19 - 15e) ^ ' 
The coefficients associated with the derivative of the static pressure are given by 

p'o = :r- {{1 + (2<?o + i^ogi) + 1} e„ , (loo) 

P'l = A {(1 + e)'^o (391 + 21/052) - 1} 00 , (101) 

P2 = ^ {(1 + e)i^i (492 + 31^093) + 1} 00 , (102) 



= A {(1 + e)'^o (593 + 41.094) - 1} So . (103) 
The coefficients associated with the derivative of viscosity are given by 

^'o = ^{29o + i^o9i)0l/\ (104) 

V ZTT 

v'o = — -9- [-17)91 + 9o {cri'^o (290 + i/09i) + t>n}] ■ (105) 

t'()9o 

The coefficients associated with the derivative of the heat conductivity are given by 

«o = —9 + So {c«fo (290 + 1^091) + b^}] , (106) 

^^090 

k'i = [a« {21^09? + 90 (91 - 21/092)} + 9o {c«!/o (391 + 2i/o92) - ■ (107) 

''o^o 

Appendix C Solution of linearized equation for the non-layering mode 

In this appendix, we solve the linearized equation for the non-layering mode {kx ^ 0) 

|^^NL^=£(%NL^ . (108) 
At first, we solve the eigenvalue problem 

= ^\llt)^^^lt) ' (109) 
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where -^j^^j) and V'k(t) = 1: 2, 3, 4) are respectively the eigenvalues and the eigenvectors 
of C{t). If we scale the wave number as k(4) = eq(t) = ^iqx,qyit)) and perturbatively solve 
Ea. lll09ll . the eigenvalues are readily found to be 

= -^^raqitf , (110) 

^qW = -^^^bQitf , (111) 

^qW = -^^r,q{tf +ieJqit) , (112) 

^qW = - , (113) 

which are respectively given by the eigenvectors 

(2) / (115) 

^q(t) ~ V2J' 2g(t)' 2 g(t) ' J 

' VQ i qx i qyjt) 
J2j'~2'q{f)'~2 q{t) ' J 



^iV. = fo,^?#,-^,0) , (114) 



(3) ^fy^ 1^ i9«(i) Po^^ , . 

^'iW V2J' 2g(t)' 2 <?(*)' jj ' ^ 



and the left eigenvectors 



im _ ^0,^,-^,0) , (118) 



V ' <?(*) ' q{t) ' 

(2) _ / 2po „ ^ 



-qW - V"~'°'°'f' ' ^^^^^ 



7,(3) _ fPo PO ^ 

V) - W' <?(*)' \(t) 'eojj 



(120) 



^(4) /J^ .9^ (121) 

where we defined q{t) = ^ g^. + qy(t)^, 

J = ^2p2/eo + ;.op^ , (122) 

and the positive constants 

_rto 2uop[^ _^ no 2pgKo , . 

- 2 ' '^^ - J2 ' - 2 + 4 + 9o J2 • ^^^^^ 

The solution of Eq. JTOSj is constructed as [19J 

= / I^G.MqW. q'Co))^.^), , (124) 

where the indexes a,/3 (= 1,2,3,4) represent the components of ■/'qj^j and we used the sum- 
mation rule for the twice appearance of /3. The Green's function is given by 



G„^(q(t),q'(0)) = EV'q{l)a'^q'(0)/3«*'^(qW'q'(0)) (125) 
with the function G'^'(q(t), q'(0)) satisfying 

{ft+"l-^^) G(^'(q(t).q'(0)) = A[;('^)G(j)(q(t),q'(0)) . (126) 



19 



Such a function G'^'(q(t), q'(0)) is found to be 

G(-')(qW,q'(0)) = (27r)25(q'(0) - q(i)) exp 



Jo 



(127) 



If we adopt y^q^^j — i V'f^^Q^ the initial condition |l9j , the components of ^^^^ 

given by 



«q(i) 
9q{t) 



.^i,(2)(t) + !|,,{3)(j)eos.W 
et „n^... 1 



1 



ii;(i){t). 



sinw(t) 
sina;(t) 



Pi)s(2)(t) + ^ij(3)(t)cosa;{t) 



where we defined 



£;(i){t) = cxp 

= exp 
£;(3)(t) = exp 



and the frequency 



QxJ 

uiit) = 



et^/l+Jd^ + In i et + ^J~^Tie^^ 



(128) 
(129) 
(130) 
(131) 

(132) 
(133) 
(134) 

(135) 



Appendix D Perturbative calculation of eigenvalue problem for the 
layering mode 



In this appendix, we perturbatively solve the eigenvalue problem of the layering mode 



£(o,fc,o)^t^^„'=-(^)(fc.o)v^^; 

where j = 1, 2, 3, 4 and the 4x4 real matrix is given by 



LO) 




C{0,ky0) 



UokyO 



_ 20.4-2 

2 "yO 



^p'okyo -(Co + Y)klo 





490 2'" 



(136) 



(137) 



\€'^ci - 2Aofc2(, -2er]okyo C2kyo e^c^ - 2Kofeyo/ 



Here, Aq ~ 0{e ). The wave number is scaled as kyo = eg and we expand C{0,kyo), a'-^'(kyo) 
and into the series of e as 



'^(O.feyo) = €Mi+e^M2 + ■ 
<f =V'^'+^^^^' + •• 



(138) 
(139) 
(140) 
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where 



Ml 



( um \ 
0-1 



-V'ol 

V 2pog o" y 



9o' 



M7 



( ° 





° \ 




2 y 










-«o + f)g2 




-2i?o(J 


C3 - 2Kog2y 



Substituting Eqs. ||138II - II139I I into Eq. ||136I I. we find the first nonzero terms 
at 0(e). From Ea. lll42ll . we find the eigenvalues 



= , 



= iJq 



The eigenvalues Eg. 11143 II are given by the eigenvectors 

</3<') = (0,1,0,0)T , 



(2) 
¥'0 



^,0,0,^ 



(3) 
¥'0 



(4) 



Bo J 

1_ l_ P0_' 

,2j'~2Jg' 2'T, 

Uo 1 « PO 

,2J' 2Jg' 2'y, 

respectively, and the corresponding left eigenvectors are given by 



^(1) - (pL 10 PO ^ 



-(2) 



-(3) 

-fio 



-(4) 



2po t'o 
- — ,0,0, — 
J J 



J 6oJJ 

^,0,i,i^') , 
J 0oJj 



(141) 

(142) 

(143) 

(144) 
(145) 

(146) 

(147) 

(148) 
(149) 
(150) 
(151) 



respectively. The eigenvectors Egs. 11144 II - II151II are orthogonal and normalized, i.e. ip^ ip^ = 
Sjk, where j,k = 1,2,3,4 and (Sj^, is the Kronecker delta. Because the critical eigenvalue is 

a real number, we are interested in the case of cr^"*^' = cr^^' = 0. However, i^q"*^^ and ipg^^ are 



degenerated, thus we rewrite Eo. 11140 II as 

(T^^^ nrirl n^^^ 

At 0(e2) of Ea. l[T36t . we find 
If we respectively multiply i^q ^ (Z = 1, 2) to Ea. lll53ll . we find 



where = (a^'^aj'')'^ and 

M2 : 



(152) 



(153) 



(154) 



(155) 
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From Eq. ||154I| . we find the eigenvalues 



7^1' - 



(156) 



which are given by the eigenvectors 



(157) 

respectively, where /i , /2 , /s and Ca are listed in Table [3] 

Because ctj^' — CTj'^' = \J — /2/2So>^^ > 0, the growth rate is given by CTj'^^ . If we expand 
(Tj^' into the series of the scaled wave number q, we find the dispersion relation 

= r2q^ + nq-^ + 0{q^) , (158) 
where the coefficients r2 and r4 are respectively given by 
_ i^O??o(2??Qpo - 7?0Pq - pQCi + 6>0PqC3) 



2!^o(PoCi - foPoCa) + 4po'?0 



(159) 



''4 = -— "7 — ^ ; — :-— ^ (ci - 2»?o)po% J + 4c3i/o6»oPo «0 (160) 

+Po{T;gj2 - 4ci!/oPoKo - %(foC3J^ + 8poKo)}] • (161) 
If we multiply ¥'o"'<^o"' (" = 3.4) to Eg. fTSSl l. we find tp^j^^ as 

^i'' = - E 4)^o"' [^<"'A^2{a(^V^^'+4'Vf }] ^C,, (0,0,1, 0)T , (162) 

n=3,4 0"l 

where 

_ 2^0)70 (1) , (aeoPoKog^ +P0C1 -6»oPoC3)po (1) 

Therefore, the eigenvector Eq.I I140I I truncated at 0(e) is given by 

vt^J) = (.„«„«;„e,)T = (^-|L4l),aW,.C^„^4l')^ . (164) 

In the same way, we calculate the left eigenvector i^^''^^ of Ea. lll36ll . Because i^q"*^' and i^q^' 
are also degenerated to the same eigenvalue tr^^' = Cj^' = 0, we write the left eigenvector as 

where / = 1, 2 and the coefficients a^'' and cij ' are determined later. At 0(e'^), we find 

-Wa^ , f-(0-(i) , -(0-(2)i {i)r-(o-(i) , -(0-{2)i 

.^i^A^i + |aiVo +4 Vo /•^2 = o-^^ jaiVo +4 Vo / ■ (166) 
If we respectively multiply i/Jq ^ (/ = 1, 2) to Eq. lll66l l. we find 

MjaC) = 4''a(0 , (167) 
where = (a^'' , )^. Then, we solve the eigenvalue problem Eq. lll67l l and find 
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Table 3 The functions fi, /2, /s, /4, /s, fe, Ca and Ca in the text. Here, Ca and Ca 
determined by the normaUzed condition of the eigenvector. 

/i = 9o{rioJ^ + 4:Uop'QKo)q'^ + irjopo + 2i'opoci - 2uo9op'oC3 

f2 = WuoO'^riop'QKoJ^q'^ + 8uoeor]oJ^{po{ci - 2ri'Q) + p'g{rio - eoC3)}q^ 

h = So(r]oJ^ - 4uoPQKo)q^ + 4»?opo - Si^oPoci + 2i^ofoPo^3 

/4 = 6'o{8poPo'^o + {2po»?o - Po'Vo)..^^}?^ + 4po(poci - eopgcs) 

fs = {uodov'oJ^ -^P()i^O + VoJ^{po+9o)}q^ +2pQ{2rio + uoci + 2poC3) 

k = /3 - (/? - /2)l/2 

Ca = -eoJ^q{{ui + 4p2)/| + 2(i.oPo - 2Bopop'o)hk + (SIjW + egPo^ + pg)/|} 



-1/2 



where /4 and Ca are listed in Table [3] If wo multiply 1,50 "^0 (" = 3:4) to Ea. |[T66ll . we find 
-(1) 

ip\ as 

^i" = - E 4) +4'Vf } ^42^"'] ^ C^i(0,0,l,0) , (169) 

n=3,4 ""l 

where 

and /s is given in TabO Therefore, the left eigenvector of Eq. lll36ll truncated at O(e^) is given 

by 

V'fcjo «9."'9.^9J - j-«2 •'^i '^"^i' 6*0 J^g ^ j ■ ' 

Let us compare the analytic form Ea. lll64ll with the numerical result. Figure |6] shows the 
components of the eigenvector {a.)uq, {h)uq, {c)wq and (d)^^ as the functions of q, where the 
open circles and solid lines represent the numerical results and analytic forms, respectively. 
From these results, Eq. ||164| I well describes the numerical results. We also confirmed that 
Ea. 1117111 well reproduces the results of the numerical calculations. 
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